// TODO: NOTE: This subroutine (except for the declaration in the header and here) has NOT been double checked.
//
// TODO: NOTE: Should update this subroutine to take in a function t (for the test statistic), so that test statistics other than the mean can be calculated.
//
// TODO: NOTE: Some of this subroutine (e.g., shuffling data, resampling, etc.) should probably make use of the resampling/bootstrapping part of stats++.


// STL
#include <algorithm>              // std::shuffle()
#include <random>                 // std::random_device, ::mt19937
#include <vector>                 // std::vector<>, ::begin(), ::end()
// stats++
#include "statsxx/statistics.hpp" // statistics::mean()


//
// DESC:
//
// =====
//
// NOTE:
//
inline double distribution::permutationTest(
                                            const std::vector<int>& X_A,
                                            const std::vector<int>& X_B,
                                            // =====
                                            const int n
                                            )
{
    // PRNG
    std::random_device rd;
    std::mt19937 g(rd());


    std::vector<int> X_AB = X_A;
    X_AB.insert(std::end(X_AB), std::begin(X_B), std::end(X_B));

/*
    for (auto& _X_A: X_A)
    {
        std::cout << "_X_A: " << _X_A << '\n';
    }
*/

    double T_obs = statistics::mean(X_A) - statistics::mean(X_B);

    double p = 0.;

    for (int i = 0; i < n; ++i)
    {
        std::shuffle(std::begin(X_AB), std::end(X_AB), g);

        std::vector<int> _X_A(std::begin(X_AB), (std::begin(X_AB) + X_A.size()));
        std::vector<int> _X_B((std::begin(X_AB) + X_A.size() + 1), std::end(X_AB));

        double T = statistics::mean(_X_A) - statistics::mean(_X_B);

//        std::cout << "T_obs: " << T_obs << " T: " << T << '\n';

        // NOTE: Comparing signed values gives a one-sided p value.
        //
        // NOTE: Comparing absolute values gives the two-sided one.
        //
        if (std::abs(T) >= std::abs(T_obs))
        {
            p += 1.;
        }
    }

    p /= static_cast<double>(n);

    return p;
}
